Constrained free-slip via recoverable Lagrange multiplier (dynamic topography)#219
Closed
lmoresi wants to merge 3 commits into
Closed
Constrained free-slip via recoverable Lagrange multiplier (dynamic topography)#219lmoresi wants to merge 3 commits into
lmoresi wants to merge 3 commits into
Conversation
…graphy Enforce u.n = g on curved boundaries with a recoverable Lagrange multiplier instead of a tuned penalty. An augmented-Lagrangian (ALG2) outer loop carries both the existing penalty BC and the multiplier; the multiplier removes the penalty's accuracy bias, so a moderate, well-conditioned augmentation gives an exact constraint in 2-3 Stokes solves. The converged multiplier is the normal traction holding the boundary -- a direct dynamic-topography estimate (h = lambda / (drho g)), recoverable as a clean MeshVariable (interior exactly zero; boundary trace correlates with -n.sigma.n at 1.0000). Purely additive: the validated 2x2 saddle-point assembly and fieldsplit config are untouched. New class behind uw.systems.Stokes_Constrained. Hands-off solve(): relative constraint tolerance (RMS(u.n-g) < rtol*RMS|u|) and viscosity-weighted augmentation (1e3*mu(x)) by default -- no user tuning. On the SolCx benchmark this gives ~2e-4 accuracy in 3 outer iterations across viscosity contrast 1 -> 1e6. Validation: tests/test_1061_constrained_freeslip_annulus.py (5 tests). Design note and production roadmap in docs/developer/design/. Underworld development team with AI support from Claude Code
Contributor
There was a problem hiding this comment.
Pull request overview
Adds a new Stokes solver variant that enforces curved-boundary normal-velocity constraints (u·n=g) using an augmented-Lagrangian outer loop with a recoverable Lagrange multiplier, enabling direct access to the converged boundary traction (dynamic-topography proxy) without penalty tuning.
Changes:
- Introduces
SNES_Stokes_Constrained(exported asuw.systems.Stokes_Constrained) withadd_constraint_bc(),multiplier(), andtopography()APIs. - Implements the augmented-Lagrangian outer iteration inside
solve()and a boundary-only multiplier update strategy via a nodal mask. - Adds annulus validation tests and a developer design note documenting formulation, trade-offs, and validation results.
Reviewed changes
Copilot reviewed 4 out of 4 changed files in this pull request and generated 3 comments.
| File | Description |
|---|---|
src/underworld3/systems/solvers.py |
Adds constrained Stokes solver, constraint BC registration, multiplier/topography accessors, and augmented-Lagrangian outer loop. |
src/underworld3/systems/__init__.py |
Exports the new solver as Stokes_Constrained. |
tests/test_1061_constrained_freeslip_annulus.py |
Adds regression/validation coverage comparing constrained vs penalty free-slip on an annulus and checks multiplier/topography behavior. |
docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md |
Design/validation documentation for the constrained free-slip multiplier approach and roadmap. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+2125
to
+2133
| point_indices = petsc_dm_find_labeled_points_local( | ||
| self.mesh.dm, | ||
| "UW_Boundaries", | ||
| getattr(self.mesh.boundaries, boundary).value, | ||
| sectionIndex=False, | ||
| ) | ||
| if point_indices is not None: | ||
| marker.data[point_indices] = 1.0 | ||
| mask = ( |
Comment on lines
+2271
to
+2277
| # Augmented-Lagrangian (ALG2) multiplier update, same r as the | ||
| # forward-problem penalty augmentation. Monotone-convergent for r>0. | ||
| # Restricted to boundary nodes so lambda stays a clean topography field. | ||
| for cbc in self._constraint_bcs: | ||
| resid = self._nodal_constraint_residual(cbc) | ||
| cbc.lam.data[cbc.mask, 0] += cbc.r_nodal[cbc.mask] * resid[cbc.mask] | ||
|
|
Comment on lines
+2224
to
+2226
| # Sample the augmentation r at the multiplier nodes once (it may be a | ||
| # spatial / viscosity-weighted expression). Used for the dual update. | ||
| for cbc in self._constraint_bcs: |
Confirmed at the PETSc level that a 3-field DM (u, p, h) with p and h grouped yields a 2-way u | [p,h] Schur split (no nested fieldsplit), de-risking the monolithic "arbitrary saddle-point constraints" direction. Captures the bounded assembly touch-points and remaining caveats for a future block-Schur PR. Underworld development team with AI support from Claude Code
- add_constraint_bc validates the boundary name and the 'UW_Boundaries' label, and rejects empty/None point sets (petsc_dm_find_labeled_points_local returns the vertex-0 sentinel np.array([0]) when the label is absent, which an `is not None` check would silently accept and corrupt the mask). - The outer loop no longer applies a final, never-solved multiplier update when the iteration cap is hit (u/p stay consistent with lambda) and now warns loudly (RuntimeWarning) on non-convergence. - Raise NotImplementedError when run on more than one MPI rank (the mask and the node-wise update are serial-only), instead of silently producing wrong results. Adds a test that add_constraint_bc rejects an unknown boundary name. Underworld development team with AI support from Claude Code
Member
Author
|
Superseded by #224. The constrained free-slip work is replaced by the in-saddle (single coupled solve) |
lmoresi
added a commit
that referenced
this pull request
Jun 8, 2026
…rsedes #219) (#224) * Add SNES_Stokes_Constrained: non-penalty free-slip + recoverable topography Enforce u.n = g on curved boundaries with a recoverable Lagrange multiplier instead of a tuned penalty. An augmented-Lagrangian (ALG2) outer loop carries both the existing penalty BC and the multiplier; the multiplier removes the penalty's accuracy bias, so a moderate, well-conditioned augmentation gives an exact constraint in 2-3 Stokes solves. The converged multiplier is the normal traction holding the boundary -- a direct dynamic-topography estimate (h = lambda / (drho g)), recoverable as a clean MeshVariable (interior exactly zero; boundary trace correlates with -n.sigma.n at 1.0000). Purely additive: the validated 2x2 saddle-point assembly and fieldsplit config are untouched. New class behind uw.systems.Stokes_Constrained. Hands-off solve(): relative constraint tolerance (RMS(u.n-g) < rtol*RMS|u|) and viscosity-weighted augmentation (1e3*mu(x)) by default -- no user tuning. On the SolCx benchmark this gives ~2e-4 accuracy in 3 outer iterations across viscosity contrast 1 -> 1e6. Validation: tests/test_1061_constrained_freeslip_annulus.py (5 tests). Design note and production roadmap in docs/developer/design/. Underworld development team with AI support from Claude Code * docs: record P'=[p,h] monolithic fieldsplit feasibility spike Confirmed at the PETSc level that a 3-field DM (u, p, h) with p and h grouped yields a 2-way u | [p,h] Schur split (no nested fieldsplit), de-risking the monolithic "arbitrary saddle-point constraints" direction. Captures the bounded assembly touch-points and remaining caveats for a future block-Schur PR. Underworld development team with AI support from Claude Code * Address Copilot review on #219: guards for label, non-convergence, MPI - add_constraint_bc validates the boundary name and the 'UW_Boundaries' label, and rejects empty/None point sets (petsc_dm_find_labeled_points_local returns the vertex-0 sentinel np.array([0]) when the label is absent, which an `is not None` check would silently accept and corrupt the mask). - The outer loop no longer applies a final, never-solved multiplier update when the iteration cap is hit (u/p stay consistent with lambda) and now warns loudly (RuntimeWarning) on non-convergence. - Raise NotImplementedError when run on more than one MPI rank (the mask and the node-wise update are serial-only), instead of silently producing wrong results. Adds a test that add_constraint_bc rejects an unknown boundary name. Underworld development team with AI support from Claude Code * Add SNES_Stokes_BlockConstrained: in-saddle multiplier for free-slip + topography Enforces u.n = g on a boundary via a Lagrange multiplier h carried INSIDE the saddle-point system (grouped u | [p,h] 2-way Schur split), in one coupled solve. The converged boundary trace h|_G = -n.sigma.n is dynamic topography directly, with no penalty/Nitsche parameter to tune. - Guarded generalisation of SNES_Stokes_SaddlePt: every multiplier block wrapped `if self._multipliers:` so the 2-field path stays bit-identical (M1 gate). - Augmented-Lagrangian conditioning of the constraint Schur (r = 1e3.mu, viscosity-weighted); combined (p,h) gauge nullspace for enclosed domains. - FMG-ready: field-index fieldsplit grouping (pc_fieldsplit_N_fields) + option mirror, so geometric multigrid works on the velocity block. - Boundary-only multiplier reduction (pin interior h) DEFAULT OFF: refined-safe via DMPlexLabelComplete closure but not yet lossless on the clone DM (degrades boundary trace ~5e-4); the lossless PetscSection-constraint path is next. - add_constraint_bc(boundary, g, normal, screening, augmentation, degree); multiplier()/topography() recovery API. - tests/test_1062: 11 cases (box + annulus; constraint / Dirichlet-match / topography / variable-viscosity to 1e6), all passing. Key finding: on hard constraints (SolCx, 4-wall free-slip + viscosity jump) the block does one mesh-independent solve where the outer-loop Uzawa needs 9-11 and explodes to 33x Dirichlet -- the two methods fail in opposite places (block: per-solve cost from the fat [p,h] Schur; outer-loop: outer-iteration count). Underworld development team with AI support from Claude Code * Lossless boundary-only multiplier reduction (default on) Constrain the interior (off-boundary) multiplier DOFs directly in the fine local PetscSection rather than via DMAddBoundary. This is both lossless (the constraint-boundary closure is correct because createClosureIndex has finalised the section) and refined/FMG-safe (no "after section creation" restriction), so the solved [p,h] Schur block carries only the boundary trace (~sqrt(ndof) DOFs). Collapses the multiplier-block DOF premium from ~3x to ~1.1x Dirichlet with no accuracy loss (box+lu reduced-vs-full velocity relL2 ~1e-9). Underworld development team with AI support from Claude Code * Promote in-saddle solver to Stokes_Constrained; remove outer-loop variant The in-saddle (single coupled solve) constrained free-slip solver is now the production path, exported as uw.systems.Stokes_Constrained. It is validated against the exact SolCx analytic solution and matches a Dirichlet free-slip reference to discretisation error, with the constraint enforced to machine precision. - rename SNES_Stokes_BlockConstrained -> SNES_Stokes_Constrained - remove the augmented-Lagrangian outer-loop solver + helper (_ConstraintBC); the Uzawa scheme is easy to reproduce in Python if wanted - single public export Stokes_Constrained (drop Stokes_BlockConstrained) - raise default augmentation_base 1e3 -> 1e4 (accuracy is r-independent; larger sits in the low-iteration plateau, well below roundoff) - tests: retarget the constrained test, add test_1062_constrained_solcx.py (free-slip vs the SolCx analytic); outer-loop annulus test removed - docs: update CONSTRAINED_FREESLIP_MULTIPLIER.md status Underworld development team with AI support from Claude Code * Address Copilot review on #224 - remove the unused _BlockConstraintBC.interior_mask and its construction (a P1 marker field + evaluate) — dead since the section-based reduction; avoidable setup cost - tests: use the public petsc_use_pressure_nullspace property (its setter resets cached nullspace state); fix the stale run-command path in test_1061 - docs/CONSTRAINED_FREESLIP_MULTIPLIER.md: rewrite the outer-loop sections to the shipped in-saddle formulation (no outer iteration; r is AL stabilisation only; augmentation_base 1e4; boundary-only reduction); fix the Files list and the option-table verdicts; mark the r-sweep table as historical outer-loop data Underworld development team with AI support from Claude Code
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds
uw.systems.Stokes_Constrained(SNES_Stokes_Constrained): enforceu·n = gon curved boundaries with a recoverable Lagrange multiplier instead of a tuned
penalty. The converged multiplier is the normal traction holding the boundary — a
direct dynamic-topography estimate
h = λ/(Δρ g).Motivated by the fragility of penalty/Nitsche free-slip on annulus/spherical
boundaries (the magnitude must be tuned against forcing strength and viscosity; too
weak gives spurious radial throughflow, too strong diverges).
How it works
An augmented-Lagrangian (ALG2) outer loop reusing the existing penalty BC: each
Stokes solve carries
t = [λ + r(u·n − g)]·n, andλ ← λ + r(u·n − g). The penaltyterm preconditions every boundary mode (fast convergence); the multiplier removes the
penalty's accuracy bias (exact constraint from a moderate, well-conditioned
r).Purely additive — the validated 2×2 saddle-point assembly and fieldsplit
configuration are untouched (honours "solver stability is paramount"). The outer loop
lives inside
solve().Hands-off
solve()(no tuning)RMS(u·n − g) < rtol·RMS|u|(defaultrtol=1e-3) — problem-independent.r = 10³·μ(x), so high-viscosityboundary regions are constrained as well as low-viscosity ones (a flat
ris not).Validation
tests/test_1061_constrained_freeslip_annulus.py, 5 tests, level_2/tier_b):RMS(u·n)=6.5e-5, velocity within 0.3% of the penalty reference, 2 outer iterations,interior λ exactly 0, boundary
corr(λ, −n·σ·n) = 0.9999.contrast 1 → 10⁶, no manual tuning.
Trade-offs and roadmap
The augmentation
ris a speed knob, not an accuracy knob (accuracy isr-independent — the AL property): too small just costs iterations, too large stiffensthe inner solve, the answer is never wrong. True-work is U-shaped with an efficient
basin
r ∈ [300, 10⁴].Deferred (see
docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md):adaptive
r₀for very coarse meshes; parallel (mask + nodal update are serial);rotation-null-space removal for both-boundaries free-slip; and the monolithic
P'=[p,h]fieldsplit ("arbitrary constraints in the saddle point") as a separatefeasibility spike.
Design note:
docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md.Underworld development team with AI support from Claude Code